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Abstract. Model reduction methods often aim at an identification of slow invariant manifolds 
in the state space of dynamical systems modeled by ordinary differential equations. We present a 
predictor corrector method for a fast solution of an optimization problem the solution of which is 
supposed to approximate points on slow invariant manifolds. The corrector method is either an 
interior point method or a generalized Gauss-Newton method. The predictor is an Euler prediction 
based on the parameter sensitivities of the optimization problem. The benefit of a step size strategy 
in the predictor corrector scheme is shown for an example. 
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1. Introduction. Models including multiple time scales arise in many applica- 
tions as e.g. combustion chemistry and biochemistry. The models under consideration 
are described by a system of ordinary differential equations (ODE). Model reduction 
methods for stiff multiple time scale ODE aim at an identification of the slow direc- 
tions in the phase space of the system. After a short relaxation phase, the dynamics 
of the ODE model are mainly determined by the slow modes. 

The slow dynamics are often described by so-called slow invariant manifolds 
(SIM). In the phase space of the dynamical system under consideration, trajectories 
bundle on those SIM being Effractors of nearby trajectories. These are hierarchically 
ordered with decreasing dimension towards stable equilibrium, and the dimension cor- 
responds to the set of slow (large) time scales. Fenichel analyzes slow manifolds in 
context of singular perturbed systems of ODE in a series of articles [11, 12, 13, 14]. 
It is shown there that for a sufficiently large separation of time scales (a sufficiently 
small singular perturbation parameter) there exists a mapping from a compact subset 
of the subspace of slow variables into the subspace of fast variables. The graph of this 
mapping is the slow manifold. 

Many model reduction methods identify or approximate such SIM. The quasi 
steady state approximation identifies the slow manifold of singular perturbed systems 
of ODE by setting the singular perturbation parameter to zero [7, 9, 35]. The dynamics 
arc only described by the differential algebraic equation that consists of the ODE 
for the slow variables and the implicit function relating the fast variables to the 
slow ones and defining the slow manifold. The ILDM method identifies a local SIM 
approximation based on an eigenvalue decomposition of the Jacobian of the right 
hand side of the ODE [29]. The CSP method identifies a local basis where fast and 
slow dynamics in the phase space are separated [23]. Many more methods for model 
reduction based on an analysis of a separation of time scales exist, see e.g. the review 
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[21]. 



In this article, wo discuss the development of an efficient application of a model 
reduction method based on optimization of trajectories to models including realistic 
detailed chemical combustion mechanisms. In Section 2, we review the optimization 
problem for model reduction. The numerical optimization method to solve this op- 
timization problem is explained in Section 3. A continuation method including the 
optimization method as corrector is introduced in Section 4. Existence of a homotopy 
path is regarded as wcill as a step size strategy for efficient progress. In Sec;tion -5, 
results of an application to example problems are shown. The article is summarized 
in Section 6 and an outlook is given. 

2. Optimization based model reduction. An approach for model reduction 

based on optimization of trajectories is raised in [24]. The idea is to identify cer- 
tain trajectories along which a criterion takes its minimum value that represents a 
mathematical characterization of slowness and attraction of nearby trajectories. This 

relaxation criterion is minimized subject to the constraints of the dynamics, the fixed 
parameters for parametrization of the SIM, and eventual conservation laws. 

2.1. Optimization problem. In [27], it is shown for two specific models that 
the method identifies asymptotically the correct SIM. Here we stick near to the for- 
mulation chosen there. We denote the ODE that describe the dynamics of the kinetic 
system with 



where D is the derivative of state vector z{t) at time t w.r.t. t. The function S* : D — >• 
M" is assumed sufficiently smooth in an open domain D G M" . Conservation relations 
arc valid for many models. Typically, these are restrictions to an initial value z{to) 
for an initial value problem with the dynamics (2.1) and represent the conservation 
of the mass of chemical elements in the system. 

A subset of state variables is usually chosen for parametrization of the SIM ap- 
proximation. In our formulation, these variables Zj{ti,), j G Ipy, which are called 
reaction progress variables or represented species, are fixed at to a certain value 
Zj', j G Xpv with |Tpv| < n. The aim of a model reduction method that allows for 
species reconstruction is the computation of the corresponding unrepresented variables 
Zj{t*), j ^ Ipv, at leading to a local representation of the SIM. 

A general formulation of the optimization problem (similar to [27]) to compute 
an approximation of a SIM is given by 



Dz(f) = S{z{t)), 



(2.1) 



min ^'(2) 



(2.2a) 



subject to 



Dz{t) 



= S{z{t)) 
= C{z{t)) 

= Zj{t*)- Zj% jeXj 

^ z{t) 



(2.2b) 
(2.2c) 
(2. 2d) 
(2.2e) 











-pv 







and 



t e [to,tf] 

t» e [to,tf] (fixed). 



(2.2f) 
(2.2g) 



A CONTINUATION METHOD IN MODEL REDUCTION 



3 



This means, a trajectory (piece) z has to be identified on an interval [to, tf] such that 
the trajectory obeys the ODE (2.2b) and necessary conservation relations (2.2c). The 
reaction progress variables are fixed to the desired values z**, j G Xpv at a point in 
time in the interval [to, if] in Eq- (2. 2d). In some applications, positivity of the 
state variables has to be demanded explicitly to guarantee physical feasibility in the 
iterative solution algorithm for problem (2.2) as it is done in (2.2e). 

2.1.1. Objective function. Various suggestions for a suitable objective func- 
tional \1/ have been made, especially in [25]. In [27], the objective function 



/■tf 

^{z) := / $ {z(t)) dt (2.3) 



'to 

with the integrand 

$(^(t)) = iiJ5(^(t)) s[zm?^ 

in the Lagrange term is used, where Jg denotes the Jacobian of S. 
This might be motivated by the estimate 



II Js Sh ^ II Js||2 \\Sh, 

where a small spectral norm || Js||2 relates to the attraction property of the SIM and 
a small 115112 relates to slowness. 

2.1.2. Local formulation. The approximation of points on the SIM computed 
as solution of optimization problem (2.2) is often used e.g. in computational fluid 
dynamics (CFD) and other applications. In this context, a large number of approx- 
imation points of the SIM for different values of the reaction progress variables is 
needed. This means, problem (2.2) has to be solved repeatedly for a range of values 
of 2:** , j e Ipv This is generally very time consuming. 

In this local formulation (in time) can be used. It is given by (cf. [20]) 

niin \\Js{z{U)) S{z{Uml (2.4a) 
^(**J)-' (**) 



subject to 



= C{z{U)) (2.4b) 
= z,(t*)-z**, jeXpv (2.4c) 
< z{U). (2.4d) 



In the case of the objective function ^ in (2.3), optimization problem (2.2) is semi- 
infinite. A discretization method such as a collocation of orthogonal polynomials or a 
shooting approach has to be applied to project it into a finite nonlinear programming 
(NLP) problem. If (2.4) is used, the optimization problem is a finite dimensional 
NLP problem, and the solution of the system dynamics (2.2b) is dispensable. This 
problem can be solved directly with standard NLP software as sequential quadratic 
programming [31] or interior point [19] methods. 
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2.2. Parametric optimization. In the optimization problem which is solved 

for model reduction purposes, there is a parameter dependence in terms of the fix- 
ation of the reaction progress variables. Therefore, we consider parametric finite 
dimensional NLP problems in the following. The standard problem is given as 

mm f{x, r) (2.5a) 

subject to 

g{x,r) = (2.5b) 

h{x,r)^0, (2.5c) 

where the objective function f : D x D ^ M., the equality constraint function g : 

D X D —i' M"^ , and the inequality constraint function h : D x D M"-^ depend on the 
parameter vector r G D with both D C M", D C M.""' open. The Lagrangian function 
can be written as 

C{x, A, fi, r) := f{x, r) - X^g{x, r) - fi^h{x, r). 

First order sensitivity results are given by the following theorem of Fiacco [16] . 
These results are used in context of real-time optimization, see e.g. [8], and nonlinear 
model predictive control, e.g. in [15, 40] , for embedding of neighboring solutions of a 
function of a parameter r. 

Theorem 2.1 (Second-order sufficient conditions [30]). Let x* be a feasible point 
of (2.5) for which the KKT conditions are satisfied with Lagrange multipliers (A*,/i*). 
Suppose further that 

v'^V^^jC{x*,X*,ii*,r)v>0 yv gC{x*,X*,h*), v^O, 

where C is the critical cone. Then x* is a strict local minimizer for (2.5). 
Proof See e.g. [17, 30]. □ 

Theorem 2.2 (Parameter sensitivity [16]). Let the functions f, g, and h in 
problem (2.5) be twice continuously differentiable in a neighborhood o/(x*,0). Let 
the second order sufficient conditions hold (see Theorem 2.1) for a local minimum of 
(2.5) at X* with r = and Lagrange multipliers X* , fi* . Furthermore, let linear inde- 
pendence constraint qualification ( LIC'Q ) he valid in (x* , 0) and strict complementary 
slackness, i.e. A*^ > i//ij(x*,0) = 0, i = 1, . . . , na. Then the following holds: 

(i) The point x* is a isolated local minimizer of (2.5) with r = 0, and the associated 
Lagrange multipliers X* , fj,* are unique. 

(ii) For r in a neighborhood of 0, there exists a unique once continuously differen- 
tiable function (a;(r), A(r), /ti(r)) satisfying the second order sufficient conditions 
for a local minimum of (2.5) such that 

(x(0),A(0),MO)) = (x*,A*,M*), 

and x{r) is a isolated local minimizer of (2.5) with Lagrange multipliers A(r) 
and n{r). 

(iii) Strict complementarity with respect to iJ.{r) and LICQ hold at x{r) for r near 0. 
Proof See [16]. □ 

Numerically, the inequality constraints (2.5c) can be treated with either an active 
set (AS) strategy or an interior point (IP) method. In an AS strategy, the AS is 
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updated in each iteration. Active constraints are considered as equality constraints, 
inactive constraints arc omitted. By contrast, the objective function is modified with 
a barrier term eliminating the inequality constraints in an IP framework. Therefore, 
we only regard eriuality constraints in the following. 

The necessary optimality (KKT) conditions (stationarity of the Lagrangian func- 
tion, feasibility, and complementarity) for problem (4.1) can shortly be written as 

K{x,X,iJ,,r) =0. (2.6) 

with e.g. for an active set method 

{V/(a;, r) - Vx g{x, r)X - Vx h{x, r)fj, 
g{x,r) 
hi{x,r), ieA{x), 

We are interested in the derivative of the solution (x* (r) , X* (r) , fi* (r)) of (2.6) 
with respect to the parameters. The implicit function theorem yields 

^{x,\,fi)K{x*,X*,iJ,*,r) Br{x,X,iJ,) = -DrK{x*,X*,iJ,*,r). 

The symbol D^/ denotes the partial derivative of / w.r.t. x. The same equation in 
matrix notation is 



[D^K BxK B^K] 

The KKT matrix 



D,.x 



-BrK. (2.7) 



[D,K DxK D^K] (2.8) 

is nonsingular if LICQ and second order sufficient optimality conditions [17, 30] are 
fulfilled. 

3. Numerical solution of the optimization problem (2.4). The optimiza- 
tion problem (2.4) is a constrained nonlinear least squares (CNLLS) problem of the 
form 

mm UF,{x)\\l (3.1a) 

subject to 

F2ix) = (3.1b) 

where the functions Fi : D ^ M"% i = 1,2, are supposed to be at least twice 
continuously differentiable in their open domain D C ffi". The problem is solved 
iteratively: Xk+i = Xk + tkdk, where dk & R"" is the increment and tk € (0, 1] the step 
length. 

In order to solve this CNLLS problem, we use a generalized Gauss-Newton (GGN) 
method as described in [6, 37], where the step direction (increment) is computed as 
solution of the linearized optimization problem 

min UF^ixk) + Ji{xk)d\\l (3.2a) 
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subject to 

F2{xk) + J2{xk)d = Q (3.2b) 

where Jj is the Jacobian matrix of Fj for i = 1,2. 

An AS strategy is used to treat inequahty constraints. For globaUzation, we 
employ a filter method [18]. The filter is implemented as in [38, 39]. The most 
important feature of a filter is the fact, that a step is accepted if it improves the value 
of the objective function or the constraint violation. The combination of a GGN 
method with a filter method for globalization of convergence is new to our knowledge. 

To prevent the Maratos effect [30], we use a second order correction (SOC) as in 
[38, 39] . The SOC increment is computed as solution of the least squares problem, cf. 
[30]: 

min \\d\\l 
d 

subject to 

F2{x + d) + J2{x)d = Q. 

If a trial point is not accepted by the filter after several reductions of the step size, 
a feasibility restoration phase is necessary. This can be done efficiently in context of a 
Gauss-Newton method. The aim of the feasibility restoration phase is the computation 
of a feasible point that is "near" to the last accepted iterate and acceptable to the 
updated filter. 

The goal to find a feasible point that is "near" to the last iterate x = Xk can be 
written in context of GGN methods as the following CNLLS problem 

min IjjS — ajjjo (3.4a) 

subject to 

F2{x)=Q (3.4b) 

Problem (3.4) can be solved iteratively with an increment dki with the new iteration 
index I as solution of a CLLS optimization problem 

min i Ijaife — a; + (i||2 (3.5a) 

d 

subject to 

F2{xk) + J2{xk)d = Q (3.5b) 

and initial value xP := x. The only difference between (3.5) and problem (3.2) is the 
objective function, matrix factorizations, e.g. of J2, can be reused. 

4. Continuation strategy. It is useful to follow the homotopy path of the 
zero of the KKT conditions in dependence of the reaction progress variables to solve 
families of optimization problems alike (2.2) for different parameters. 

In the following, we consider the finite dimensional parametric optimization prob- 
lem 



min fix) 



(4.1a) 
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subject to 

= g{x) (4.1b) 
= (,,;-| — r\ i = l,...,nr (4-lc) 

where the functions / : £> — ^ R and g : D ^ K"^ are C^{D) in the open domain 
£) C M". The fixed values of the reaction progress variables in (4.1c) are denoted by 
the parameter vector r e M"'', r' = i = 1, . . . , rir, <n — n2 with the notation 
of Equation (2. 2d), where j : {1, . . . , n,.} Tpv, i i-> j(i) is a bijective map to the 
index set Xpv C {1, . . . , n} of the reaction progress variables. 

4.1. Parameter sensitivities in the GGN method. If Newton's method 
together with an active set method is used to find a candidate solution of (4.1), 
i.e. the computation of a root of JC, the KKT matrix Ti(^x \^^-^K{xk,\k^^k) has to 
be computed in every iteration k. This is different if a generalized Gauss-Newton 
method is employed. 

We return to the notation used in Section 3: The objective function is written 
as f{x) = \{Fi{x)\{^ with a sufficiently smooth function Fi : D C K" -S> M"i. The 
constraints are given as F2 : D cM" ^ M"^, n2 = n2 + + |>l(a;)| with 

(9{x) 

F2{x) = I -r\ i = l,...,nr 
[^Xi, i € A{x) 

and the active set A{x) in x. The Jacobian matrices of Fi and F2 are denoted as Ji 
and J2, respectively. 

In every iteration, the equation system 



'JIJl 


•J 2 




■ d ' 













-A 







(4.2) 



is solved in the GGN method (where the argument Xk is omitted and the vector A is 
used for all equality and active inequality constraints); that is the KKT system of the 
CLLS problem (3.2), see Section 3. By contrast. 







'Ax 






J2 




AA 




. ^2 . 



(4.3) 



is solved if Newton's method is applied to find a KKT point of the original CNLLS 
problem (3.1). The derivative of the Lagrangian function of the CNLLS problem 
is Vx^^ = Fi{x)'^ Ji{x) — X^J2{x). The difference in the KKT matrices in Equa- 
tions (4.3) and (4.2) is the difference in the Hessians of the different Lagrangian 
functions: In the GGN method, J^Ji is used; whereas Newton's method applied to 
the KKT conditions of the original CNLLS problem (3.1) exploits 

"2 

= V^4l|Fi(a;)||2 - ^ A,V,, F^(a;) 

_ (4.4) 

"2 ^ ' 

= J^Ji + (D,jT) Fi -Y,X,VxxF^{x), 
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where F2 is the i-th component of F2. 

In our appHcation in chemical kinetics, the constraints F2 arc given via conser- 
vation relations. In this case, the second derivative S/xxFi oi typically is zero. 
This might not be true, if an entry in F2 is a nonlinear equation in x representing 
an energy conservation law, see also [26]. This means, the effort for computing Cxx 
instead of Ji J\ is mainly the effort to compute the expression {PxJi) Fi. This can 
be evaluated using automatic differentiation [22] as a directional derivative of Ji with 
respect to direction F\. 

4.2. Euler predictor. If the approximation of points on the SIM is used in 
simulations in computational fluid dynamics, approximations of the tangent vectors 
of the SIM in these points are needed, too, see e.g. [34]. These are given via the 
parameter sensitivities D^iX* = '^^ at the solution x* and z* of the optimization 

problems (4.1) and (2.2), respectively. As these derivatives are available, we use the 
Euler prediction in a predictor corrector scheme for initialization of the optimization 
algorithm (corrector) to solve neighboring problems. In our case, this is the compu- 
tation of a solution with the optimization algorithm for various parameter values of 
the reaction progress variables r. 

The prediction can be used in a homotopy method with a step size strategy. An 
effective step size strategy is published in [10] and extensively discussed and modified 
in [3, 4, 5]. We use the method as discussed in [4]. The aim of the step length strategy 
of den Heijer and Rheinboldt [10] is to achieve a desired number of iterations for the 
corrector step. The strategy allows for the computation of a step length based on the 
contraction rate of the latest corrector iterations and an error model for the corrector 
such that the desired number of iterations is achieved. 

We define a curve c as a mapping from the parameter space to the space of the 
primal and dual variables of (4.1) 

For each value of the parameter vector r, c(r) is the solution {x*{r), A*(r)) of (4.1). 
We denote the Euler predictor step 

c°+i{hi) = c, + h,{r,+i-r, f^.cM, i = 0, . . . (4.5) 

The prediction c° is used as initialization of the corrector to compute the solution of 
(4.1). The corrector iterations are denoted 

4+\h,) = C{c>iK)), j=0,... 

with the corrector C. It is assumed that the corrector iterations converge to the 
solution 

c~(/ii) := lim ciihi) e K-^{0). 

The sophisticated aspect in the work of den Heijer and Rheinboldt [10] is the 
error model </>. This error model estimates the error of the iterates 



ej{hi) = \\cr{hi)-4ihi)h 
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Table 5.1 

Simplified mechanism as used in [33]. Collision efficiencies M; an = 1-0, = 2.5, aoH = 
1.0, ao = 1.0, aHsO = 12.0, = 1-0; [33]. 



Reaction 






A 1 (cm, mol, s) 


h 


1 kJmori 


O + H2 




H + OH 


5.08 X 10°'' 


2.7 


26.317 


H2 + OH 




H2O + H 


2.16 X 10°^ 


1.5 


14.351 


+ H2O 




2 OH 


2.97 X 10°^ 


2.0 


56.066 


H2 + M 




2H + M 


4.58 X lO^^ 


-1.4 


436.726 


+ H + M 




OH + M 


4.71 X 10^^ 


-1.0 


0.000 


H + OH + M 




H2O + M 


3.80 X 10^2 


-2.0 


0.000 



independently of h via an expression of the form 

ej+i(/ii) < <P{ej{hi)). 

The formula for the error model depends on the contraction rate of the corrector, see 
[4, p. 53fF.], where we use the error model for the quadratically convergent Newton's 
method and for the Unearly convergent GGN method. 

Linear step. Optimization problem (4.1) has to be solved many times for different 
values of the parameter r . Especially if the approximation of points on a SIM is needed 
in situ, e.g. in a CFD simulation, it is necessary to compute the solution of neighboring 
optimization problems fast. 

If the values r'^''^ for the reaction progress variables, for which a SIM approxima- 
tion is needed, are near to the values r* 

||r"^"-r*||2<etoi (4.6) 

for which the optimization problem is already solved, it can save computing time to 
use z^™ as approximation of the SIM. This is defined (in analogy to (4.5)) as 

:= z*it,) + (r"™ - r*f ^^^(r*) (4.7) 

where the notation is in accordance with the notation in (2.2), r' = i = 1, ■ ■ ■ , n^, 
rir = |2pv|) j is the bijcction, and r* is the z-th component of r. 

5. Numerical results. For numerical validation of our method, we use a test 
model for model reduction purposes. The reaction mechanism is given in Table 5.1. 
We use thermodynamical data in form of coefficients of NASA polynomials we received 
from J. M. Powers and A. N. Al-Khateeb, which they use in [1, 2]. 

The mechanism is published originally in [28]. The simplified version shown in 
Table 5.1 is used by Ren et al. in [33]. The mechanism consists of five reactive species 
and inert nitrogen, where in comparison to a full hydrogen combustion mechanism the 
species O2, HO2, and H2O2 are removed. The species are involved in six Arrhenius 
type reactions, where three combination/decomposition reactions require a third body 
for an effective collision. The reaction system is considered under isothermal and 
isobaric conditions at a temperature of T = 3000 K and a pressure of p — 101 325 Pa. 
Hence, the state of the system is sufficiently described by the specific moles of the 
chemical species Zi, i = 1, . . . ,nspec in molkg"^. 
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Conservation relations for the elemental mass in this model are given in terms of 
amount of substance as [33] 

riH + 2 + noH + 2 riu^o = 1-25 x 10"^ mol 

«oH + no +nH20 = 4.15 X 10~^mol (5.1) 
2nN2 = 6.64 X 10"^ mol. 

The total mass in the system can be computed with the values in Equation (5.1) and 
has a value of m = 1.01 x 10~^ kg. 

5.1. One-dimensional SIM approximation. Results of an approximation of 
a one-dimensional SIM are shown in Figure 5.1 and 5.2. We use zhsO as reaction 
progress variable and vary its value. The values of all remaining species at the solution 
of the different optimization problems are plotted versus zhjO in form of trajectories 
through the solution points z*{t») which are shown as x marks. The computation is 
done with the GGN method as described in Section 3. 

We use values near the equilibrium state as initial values in the optimization 
algorithm as we assume this point near a slow manifold. We set the values of the 
reaction progress variable near its initial value; see Table 5.2. This allows for a fast 
computation of a solution of the first optimization problem with z^^^q = 3molkg~^. 

Table 5.2 

Initial value (unsealed) for the algorithm and a solution of the optimization problem (2.2) as 
solved for the results depicted in Figure 5.1 to reduce the simplified hydrogen combustion model with 
^HsO = 3molkg-i. 



Variable 


Initial value z°{t^,) 


Numerical solution z*{t^) 


zo 


0.345 46441 


0.345 637 63 


■2H2 


2.027973 2 


2.0281615 


zh 


1.519 563 9 


1.519 3606 


ZQH 


0.764 549 59 


0.764 376 37 




3.000 000 


3.000 0000 


^N2 


32.905 130 


32.905 130 



All computations are done on an Intel® Core'^'^ i5-2410M CPU with 2.30 GHz, 
operating system openSUSE 12.2 (x86_64) including the Linux kernel 3.4.11 and GCC 
4.7.1. We do not use a step size strategy in the predictor corrector scheme. To compute 
the 17 points shown in Figure 5.1, 74 iterations are necessary, which take in sum 0.02 s. 

Figure 5.2 is presented to illustrate the linear step approximation. For the com- 
putation of the three solutions (as etoi = 1-1 is chosen) of the optimization problem 
with different values for ZfjaOj 15 iterations are needed (in sum), which is slightly more 
effort per optimization problem than in the case, where the results arc illustrated in 
Figure 5.1, where the distance between the different values for z^' q is only 0.25. It 
can be seen that the result of the linear approximation can lead to large deviations 
from a smooth, invariant manifold. 

Warm, start with step size strategy. We want to demonstrate that a step size 
strategy for the warm start of the algorithm to solve neighboring optimization prob- 
lems can have a benefit. So wc consider the optimization problem (2.2) with = tf, 
ti — to = 10~^s to be solved for the nominal parameter z^ q = 3molkg~^ with the 
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2.5 



1 2 3 

Z H20 



1.5 



0.5 




1 2 3 

'Z H20 




1 2 3 

'Z H20 



1 2 3 

2; H20 



Fig. 5.1. Numerical solutions of the optimization problem (2.4) to approximate a one- 
dimensional manifold in the state space of the six species model for hydrogen combustion. The 
full dot represents the equilibrium. The x marks depict solutions z*(t») of (2.4) for different values 
of the reaction progress variables 2H2O 'J* ■ Trajectories through these points ( curves in the figures) 
are also shown. 



shooting approach and Ipopt [39]. As neighboring problem we consider the parame- 
ter %*20 = O.Smolkg^^. Such a large change in the parameter can occur e.g. if the 
presented method for model reduction is used in situ in a CFD simulation and grid 
refinements are performed in different regions of the spatial domain. 

We solve the optimization problem first with 2^*20 ^ Smolkg^^ and second with 
Q = 0.5molkg~^ with a full step method in the predictor corrector scheme and 
apply the Euler prediction. 

The computations are done on an Intel® Core'^'^ i5-2410M CPU with 2.30 GHz, 
operating system opcnSUSE 12.2 (x86_64) including the Linux kernel 3.4.11 and GCC 
4.7.1. Six iterations in Ipopt [39] have to be preformed to solve the nominal opti- 
mization problem and nine iterations to solve the second optimization problem. In 
sum, 15 iterations are necessary which take 0.82 s in sum. 
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Fig. 5.2. Illustration of solutions of the same problem as solved for the results shown in 
Figure 5.1, but with a linear step tolerance o/etol = 1-1> see (4-6). 



If we initialize the step size strategy, see Section 4.2, with initial step size hmit = 
and the desired number of iterations k — 10, we need one intermediate step in the 
predictor corrector scheme described in Section 4.2. This is to solve the optimization 
problem with the parameter z^I^q = 2 molkg"^. In this case we need in sum 6+5+8 = 
19 iterations that take only 0.70 s. The difference in the computation time arises as 
the point for initialization of the algorithm to solve the second optimization problem 
in the full step method is not near the solution such that the KKT matrix is ill- 
conditioned. The inertia correction in Ipopt, see [39], is activated, and 19 line search 
iterations are performed in sum. In case of an activated step size strategy, the initial 
value for the algorithm to solve the optimization problem in case of a warm start is 
near the solution of the optimization problem such that no inertia correction and also 
19 line search iterations (one per Newton iteration) are necessary in the presented 
example. 

For an efficient tracking of the SIM for different values of the reaction progress 
variables, it is important to stay close to the solution in neighboring optimization 
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problems. Therefore, a step size strategy is beneficial. 

5.2. Two-dimensional SIM approximation. We choose zi = zh^ and Z2 = 

zq as reaction progress variables. Numerical solutions of (2.2) for the reduction of 
the simplified model for hydrogen combustion are shown in Figure 5.3. 




Fig. 5.3. Visualization of numerical solutions of the optimization problem (2.2) with = if, 
t{ — to = 10"'' s to reduce the simplified hydrogen combustion model with the reaction progress 
variables z\ and z^. The solution points z*{tt) with 2:3 = zhjO /f"" different values of (z'^' jZ^') 
are shown as a mesh. The results shown in Figure 5.1 are again plotted as x marks as well as the 
equilibrium is shown as full dot. 

Local solutions of the optimization problem (2.2) approximate different two-di- 
mensional SIM for different values of the reaction progress variables. The two two- 
dimensional SIM come close to each other. Such occurrences can lead to severe numer- 
ical problems in the optimization algorithms as well as in the continuation algorithms 
as there are regions where the KKT matrix might be singular at least in the range of 
machine precision. 

5.2.1. Optimization landscapes. We want to further illustrate this situation 
via optimization landscapes, i.e. graphical representations of the objective function 
versus the (free) optimization variables. 

One reaction progress variable. In Figure 5.4, the value of (the objective function 
of the optimization problem (2.4)) $ = ||Js(2;) ^(z)!!! is plotted versus the two 
degrees of freedom in (2.4) represented by zi and Z2 for one exemplary value of the 
one reaction progress variable Zg* = 2:^*20 ^ Imolkg"^ for the simplified hydrogen 
combustion model. 

The $-axis has logarithmic scale. A single local solution of the optimization 
problem (2.4) for the reduction of the simplified hydrogen combustion model can be 
seen. 

Two reaction progress variables. To compute an optimization landscape in case 
of two reaction progress variables, we fix — = O.Bmolkg"^. We regard the 
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Fig. 5.4. Visualization of ^ = \\Js{z) S{z)\\2 in dependence of z\ and Z2 to illustrate the 
solution of the optimization problem (2.4) with z^* = z^'^q ~ Imolkg^-'- for the reduction of the 
simplified hydrogen combustion model. 

value of $ = \\Js{z) S'(z)||2 (the objective function of the optimization problem (2.4)) 
in dependence of Z3 for fixed zi. 

The results are shown in Figure 5.5. It can be seen that there are two distinct 
local minima of $ = \\Js{z) S{z)\\2 for a fixed value of e.g. zi = 2molkg~^: There 
is no unique local solution of the optimization problem (2.4) for the reduction of 
the simplified hydrogen combustion model with the reaction progress variables = 
Zq = 0.3molkg~^ and 2** = z^*^ = 2molkg^^. One solution is near the value 
Z3 = z^^Q — 3.1020molkg^^ and the other one near Zg — 5.2977 molkg"^. 

In this case the predictor corrector scheme helps to follow the desired local optimal 
solution in dependence of the reaction progress variables and not to switch to another 
curve c of local solutions. This is only possible as long as the KKT matrix is not "too 
ill-conditioned" . 

5.2.2. Performance test. We use the specific moles of H2O and II2 for the 

parametrization of a two-dimensional SIM approximation in a performance test. 

We consider a test situation of a two-dimensional grid of 108 points (zjj'aO' "^112) ^ 
[0.001, 0.5, 1, 1.5, . . . , 5, 5.5] x [0.001, 0.5, 1, . . . , 3.5, 4], where points which violate mass 
conservation in combination with the positivity constraints are ignored such that 80 
points remain. 

Solutions of the optimization problem (2.4) computed with the GGN method with 
Euler prediction for initialization of the algorithm to solve neighboring problems are 
shown in Figure 5.6 for illustration. 

In Table 5.3, we compare the performance of the algorithm using the different 
implemented solution methods. We use the generalized Gauss-Newton method as 
described in Section 3 for the solution of (2.4). We solve the same problem with the 
interior point algorithm Ipopt [39]. Third, we apply a shooting approach for the semi- 
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Fig. 5.5. Visualization of ^ = || J5(2) ^(z)!!! for the simplified hydrogen combustion model 
in dependence of zi and 23 for a fixed value of Z2 = 0.3molkg~^. The scale for $ is logarithmic. 
Consider the value of ^ in dependence of z-j for a fixed value of zi, e.g. for z\ = 2molkg~^. 

infinite optimization problem (2.2) with = t{. The NLP problem (after solving the 
ODE with the BDF integrator developed by D. Skanda in [36]) is solved with Ipopt 
[39]. As a fourth alternative, we use a simple Gauss-Radau collocation with linear 
polynomials (backward Euler) for (2.2) with ~ t{. The resulting high-dimensional 
NLP problem is also solved with Ipopt [39] in the latter case. For the optimization 
problem (2.2), an integration horizon of <t — = 10~^ s is used. 

Table 5.3 

Comparison of the performance of the various algorithms for the reduction of the simplified 
model for hydrogen combustion with two reaction progress variables and full step predictor corrector 
method for the initialization of neighboring problems. 



Method 


Prediction 


^ Iter, w/o fail 


Time 


Fail 


Time w/o fail 


GGN 


Constant 


806 


0.43 s 


11 


0.18s 




Euler 


731 


0.35 s 


11 


0.17s 


IP for (2.4) 


Constant 


670 


1.09 s 





1.09s 




Euler 


591 


0.85 s 





0.85 s 


Shooting 


Constant 


335 


44.15 s 





44.15 s 




Euler 


258 


39.36 s 





39.36 s 


Collocation 


Constant 


315 


95.61s 





95.61s 




Euler 


234 


51.13s 





51.13s 



An initialization of neighboring problems is done without step size control, as we 
evaluate the benefit of the Euler prediction here. The computations are done with 
the same computer configuration as in the example for a warm start with step size 
control in Section 5.1. 
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Fig. 5.6. Illustration of an approximation of a two- dimensional SIM in the phase space of 
the simplified model for hydrogen combustion. The optimization problem (2.4) is solved with the 
GGN method. The absolute tolerance for all di is 10"^''. The relative tolerance is 10~^. The 
values of the optimization variables at the solution are plotted versus the given values of the reaction 
progress variables shown as x marks. The trajectories started with the solution z*{tt) as initial 
values converge toward the eguilibrium shown as full dot. 



It can be seen that the gain in the computation time achieved by the Euler 
prediction in case of the collocation method, where we use 100 collocation points, is 
the largest with about 46.5%. This is reasonable because of the strong dependence 
of the collocation method on the initialization on all collocation points in the time 
interval [io, ^f]- 

The benefit of the Euler prediction is about 10.8% if the shooting approach is 
used to solve (2.2). In case of the GGN method, the benefit is only 5.6% if we do 
not regard failures. The eleven failures only occur in the region, where z^* q is larger 
than the equilibrium value z'^ q. These can not be overcome neither with the step 
size control for the continuation method nor with a larger tolerance for convergence 
in the GGN method. It can be seen that the algorithm for the solution of (2.4) with 
the GGN method is much faster than the algorithm for the solution of (2.2). The 
computation of one approximation of a point on the SIM (without regarding failures) 



A CONTINUATION METHOD IN MODEL REDUCTION 



17 



takes about 0.17s/(80 — 11) « 2.5 ms. This is pretty fast and migiit be used in an 
online {in situ) SIM computation during CFD simulations. 

6. Conclusion. In this article, we present a strategy for an efficient solution of 
parametric optimization problems that arise for a method for model reduction that 
is raised in [24, 25, 32]. 

Two methods are tested to solve the semi-infinite optimization problem (2.2) 
for model reduction: collocation and shooting approaches. The resulting nonlinear 
programming problem is solved with a state-of-the-art open source interior point 
algorithm [39]. It turns out that all variants for the solution of the semi-infinite 
optimization problem arc too slow for an in sUm application of the model reduction 
method in e.g. a computational fluid dynamics simulation. 

A finite optimization problem in form of a constrained nonlinear least squares 
problem can be formulated instead. This can be solved efficiently with a generalized 
Gauss-Newton method. A filter approach is used for globalization of convergence. 
Second order correction iterations prevent the Maratos effect. The problem that has 
to be solved in the feasibility restoration phase of the filter algorithm can also be 
formulated as a least squares problem such that matrix factorizations can be reused. 

Parameter sensitivities of the optimization problem are used in a homotopy 
method for the solution of neighboring problems. The reaction progress variables 
for parametrization of the slow manifold are considered as parameters in the opti- 
mization problem. The problem has to be solved several times for different values 
of the parameters. An Euler prediction of the solution of neighboring problems is 
employed based on the parameter sensitivities, which are computed to obtain an ap- 
proximation of the tangent space of the slow manifold. This linear prediction can be 
used directly within a certain tolerance as an approximation of a point on the slow 
manifold. It can also be used in combination with a step size strategy. The step size is 
computed according to [10] in dependence of the contraction of the corrector method 
- the solution method for the nonlinear programming problem - and the iterations 
needed to solve a previous optimization problem. 

We study the presented methods for a simplified model for hydrogen combustion. 
We find that a solution of the optimization problem for approximation of points on a 
slow manifold can be computed in short time for the presented example. Furthermore 
it can be seen that there might exist several distinct local solutions of the optimization 
problem. In such cases, the predictor corrector scheme helps to follow one local 
minimum in dependence of the chosen parametrization of the slow manifold. 

Such a solution strategy could also be used for variations in the model parameters 
as e.g. the mixture fraction, the internal energy, and others. 
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